function y = func_delta(p, m, CONSTS)
    
    k0    = CONSTS.k0;
    a     = CONSTS.a;
    eta   = CONSTS.eta;
    eps_a = CONSTS.eps_a;

    s      = func_s(p, CONSTS);
    n1     = func_n1_from_p(p, CONSTS);
    n2     = func_n2_from_p(p, CONSTS);
    K_m    = func_K_m(p, m, CONSTS);
    J_m1   = func_J_m1(p, m, CONSTS);
    J_m2   = func_J_m2(p, m, CONSTS);
    J_m1_1 = func_J_m1_1(p, m, CONSTS);
    J_m2_1 = func_J_m2_1(p, m, CONSTS);


    y = (n2-n1).*((p.^2)./eps_a).*((m^2)./((k0*a*s).^4)) - (n2-n1).*(K_m).^2 + ...
        n1.*(eta/eps_a).*J_m1_1.*J_m2_1 - n2.*(eta/eps_a).*J_m1.*J_m2 - ...
        n2.*J_m1.*K_m + (eta/eps_a)*n1.*J_m1_1.*K_m + n1.*J_m2_1.*K_m - ...
        (eta/eps_a)*n2.*J_m2.*K_m + p.*(eta/eps_a).*(m./((k0*a*s).^2)).* ...
        (-J_m2_1 + J_m1 + J_m1_1 - J_m2);

end

function y = func_J_m1(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q1 = func_q1_from_p(p, CONSTS);
    alfa_1 = func_alfa_1(p, CONSTS);
    bessel_m_q1 = besselmx(double('J'), m, k0*a*q1, 0);
    bessel_mp1_q1 = besselmx(double('J'), m+1, k0*a*q1, 0);

    y = bessel_mp1_q1./((k0*a*q1).*bessel_m_q1) + ...
           (alfa_1.*m./((k0*a*q1).^2));
end

function y = func_J_m2(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q2 = func_q2_from_p(p, CONSTS);
    beta_2 = func_beta_2(p, CONSTS);
    bessel_m_q2 = besselmx(double('J'), m, k0*a*q2, 0);
    bessel_mp1_q2 = besselmx(double('J'), m+1, k0*a*q2, 0);

    y = bessel_mp1_q2./((k0*a*q2).*bessel_m_q2) - ...
           (beta_2.*m./((k0*a*q2).^2));
end

function y = func_J_m1_1(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q1 = func_q1_from_p(p, CONSTS);
    beta_1 = func_beta_1(p, CONSTS);
    bessel_m_q1 = besselmx(double('J'), m, k0*a*q1, 0);
    bessel_mp1_q1 = besselmx(double('J'), m+1, k0*a*q1, 0);

    y = bessel_mp1_q1./((k0*a*q1).*bessel_m_q1) - ...
           (beta_1.*m./((k0*a*q1).^2));
end

function y = func_J_m2_1(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q2 = func_q2_from_p(p, CONSTS);
    alfa_2 = func_alfa_2(p, CONSTS);
    bessel_m_q2 = besselmx(double('J'), m, k0*a*q2, 0);
    bessel_mp1_q2 = besselmx(double('J'), m+1, k0*a*q2, 0);

    y = bessel_mp1_q2./((k0*a*q2).*bessel_m_q2) + ...
           (alfa_2.*m./((k0*a*q2).^2));
end

function y = func_K_m(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    s = func_s(p, CONSTS);

    y = besselk(m+1, k0*a*s, 0)./((k0*a*s).*besselk(m, k0*a*s, 0)) - ...
           (m./((k0*a*s).^2));
end
